#include "AnaDef.h"

#define EPSILON 1e-10

TH1D * getNewNSigma(const TH1D * hsignal)
{
  const TAxis * ax = hsignal->GetXaxis();
  TH1D * hnsigma = new TH1D(Form("%snsigma", hsignal->GetName()),"", ax->GetNbins()+1, ax->GetXmin()-ax->GetBinWidth(1)*0.5, ax->GetXmax()+ax->GetBinWidth(1)*0.5);

  const Double_t totn = hsignal->Integral(0, hsignal->GetNbinsX()+1);
  Double_t hmin = 1e10;

  //this is the sigma in the definition of TMath::Erfc
  const Double_t sigma = 1/TMath::Sqrt(2);

  if(totn>EPSILON){
    for(Int_t ii=1; ii<=hnsigma->GetNbinsX(); ii++){
      const Double_t frac = hsignal->Integral(ii, hsignal->GetNbinsX()+1)/totn; 

      //avoid ==0, ==1 for numerical constraint in TMath::ErfcInverse
      if(TMath::Abs(frac-1)<EPSILON || frac<EPSILON){
        hnsigma->SetBinContent(ii, -999);
        continue;
      }

      const Double_t arg = frac * 2;
      //http://root.cern.ch/root/html/src/TMath.cxx.html#m5Z1cC
      //printf("test ii %d totn %e frac %e %e arg %e\n", ii, totn, hsignal->Integral(ii, hsignal->GetNbinsX()+1), frac, arg);
      const Double_t xsigma = TMath::ErfcInverse(arg)/sigma;
      hnsigma->SetBinContent(ii, xsigma);      
      hnsigma->SetBinError(ii,0.01);

      if(xsigma<hmin){
        hmin = xsigma;
      }
    }
  }
  
  if(hmin>0)
    hmin*=0.9;
  else
    hmin*=1.1;

  hnsigma->SetMinimum(hmin);
  return hnsigma;
}

void demo(const Int_t id=0, const Int_t ncount = 10000)
{
  TH1D * hsignal=new TH1D("hsignal","",100,0,10);

  TRandom3 ran(111);
  TString tit;
  for(Int_t ii=0; ii<ncount; ii++){
    Double_t val = -999;
    switch(id){
    case 0:
      val = ran.Gaus(5,1);
      tit="Gaus";
      break;
    case 1:
      val = ran.Landau(5,0.1);
      tit="Landau";
      break;
    case 2:
      val = ran.PoissonD(5);
      tit="PoissonD";
      break;
    case 3:
      val = ran.Exp(0.5);
      tit="Exp";
      break;
    case 4:
      val = ran.Rndm()*4+3;
      tit="Rndm";
      break;
    case 5:
      val = ran.BreitWigner(5,1);
      tit="BreitWigner";
      break;
    default:
      printf("bad id %d\n", id); exit(1);
    }

    hsignal->Fill(val);
  }
  hsignal->SetTitle(tit);

  gStyle->SetOptStat(0);
  TCanvas *c1=new TCanvas("c1","",400,600);
  c1->Divide(1,2);
  c1->cd(1);

  hsignal->SetYTitle("Counts");
  hsignal->Draw();
  if(id==0){
    hsignal->Fit("gaus");
  }

  c1->cd(2);
  gPad->SetGrid();
  TH1D *hnsigma = getNewNSigma(hsignal);
  hnsigma->SetYTitle("N#sigma");
  hnsigma->SetTitle(tit);
  hnsigma->Draw("hist");
  if(id==0){
    hnsigma->Fit("pol1");
  }

  c1->Print(Form("ns%d.eps", id));
}

int main(int argc, char *argv[])
{
  for(Int_t ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }
  if(argc!=3){
    printf("bad argc!\n"); return 1;
  }

  demo(atoi(argv[1]), atoi(argv[2]));

  return 0;
}
